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Abstract 

The role played by the composite analogue of the log likelihood ratio in hy- 
pothesis testing and in setting confidence regions is not as prominent as it is in the 
canonical likelihood setting, since its asymptotic distribution depends on the un- 
known parameter. Approximate pivots based on the composite log likelihood ratio 
can be derived by using asymptotic arguments. However, the actual distribution 
of such pivots may differ considerably from the asymptotic reference, leading to 
tests and confidence regions whose levels are distant from the nominal ones. The 
use of bootstrap rather than asymptotic distributions in the composite likelihood 
framework is explored. Prepivoted tests and confidence sets based on a suitable 
statistic turn out to be accurate and computationally appealing inferential tools. 
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1 Introduction 



1 . 1 Overview 

When dealing with complex models, canonical likelihood inference may encounter some 
theoretical and computational difficulties. For instance, in models with complicated tem- 
poral and/or spatial dependence structures, a likelihood function based on the joint dis- 
tribution of the observable data might even be unavailable. On the other hand, the 
specification of the joint distribution can be straightforward, but the evaluation of the 
likelihood function might lead to computational burden. In order to cope with these 
difficulties both in model specification and in computation, the use of composite likeli- 
hood functions may prove useful, as advocated by several authors both in the frequentist 
domain (see, e.g., Varin et al, 2011) and, more recently, in the Bayesian setting (see. 



e.g., Pauli et al.. 2011). Composite likelihoods have shown a great impact also in prac- 



tical applications. Some examples are spatial processes (Varin et al., 2005), multivariate 



extremes (Padoan et al.. 2010), and longitudinal models (Fieuws and Verbeke, 2006). 



In spite of the high flexibility and multiplicity of applications of composite likelihood 
functions, some concerns arise about the accuracy of the derived inferential procedures 
when testing and constructing confidence sets for a multidimensional parameter, as the 
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use of the composite log likelihood ratio is not as straightforward as it is in the canonical 
likelihood setting. In fact, its asymptotic distribution is non-standard and depends on 
unknown coefficients that need to be estimated from a matrix related to the Godambe 



information (Kent, 1982). Tests and confidence sets can also be defined by considering 



the usual Wald and score statistics as well as on suitable modifications of the composite 



log likelihood ratio (Geys et al, 1999 Chandler and Bate 2007; Pace et a/. 2011). 



Nevertheless, the evaluation of the aforementioned statistics also requires the computation 
of the Godambe information. Therefore, the accuracy of composite likelihood inference 
relies upon the Godambe information matrix that, as a matter of facts, regulates the rate 
of convergence of the sampling distribution of statistics to their asymptotic references. 

Inference based on asymptotic approximations can be improved by resorting to boot- 
strap techniques. Aerts and Claeskens ( 1999 ) propose the use of parametric bootstrap 
to approximate the distribution of general pseudo-log likelihood ratios. However, this 
approach can be computationally intensive, and as a major drawback, it requires the 
specification of the joint distribution of the data. Also the semiparametric bootstrap 
could be considered, but its application is limited to a narrow range of applications 
(Aerts and Claeskens, 2001). 

The aim of this work is to motivate the use of nonparametric bootstrap in the com- 
posite likelihood framework. Stemming from the original formulation of prepivoting in- 
troduced by Beran (1987, 1988) and refined by Lee and Young (2003), bootstrap theory 
developed in standard settings is conveyed to models involving highly structured depen- 
dencies. Prepivoting has been proven to be a general and effective approach alternative 
to the use of asymptotic refinements that allows to reduce the error level of tests and 
confidence sets. However, it has been largely neglected because its application usually 
requires a computationally expensive Monte Carlo simulation. It is shown how prepivot- 
ing a suitable statistic, namely the unstudentized quadratic form of the composite score 
statistic, aids at circumventing the computational difficulties and at the same time yields 
accurate inferential procedures. 

In the remainder of this section composite likelihood functions are reviewed, espe- 
cially with reference to marginal pairwise likelihood functions, and a general formulation 
of prepivoting is outlined and contextualized in the pairwise likelihood framework. A 
description of the proposed prepivoting approach is presented in Section [2] and its finite 
sample properties are assessed via Monte Carlo simulation in Section [3j Finally, a brief 
discussion is given in Section |4j 



1.2 Marginal pairwise likelihoods 
1.2.1 Definitions and notation 

In the following, denote with y = {yi, . . . , y„) a sample of independent realizations of the 
random vector y G R'^ supposed to have probability distribution Fg and density function 
f{-]0) depending on a multidimensional parameter 6* G O C Rp. Let i{6) = log f{y] 6) 
be the log likelihood function and w{9) = 2[i{9) — i{9)] be the log likelihood ratio, with 
6 the maximum likelihood estimate. 

Consider a set of marginal measurable events {Sj. E y, r = 1, . . . ,m} defined on the 
sample space 3^ and let fr{yi',9) = f{yi G Sr;0), i = 1, . . . , n, be the likelihood contribu- 
tion generated from f{yi; 9) by considering the set £r- The composite likelihood function 



2 



is defined as the product of sub-likelihoods 



=1 r=l 



where Ur are non-negative weights. 

The marginal pairwise likelihood function is a subclass of composite likelihoods ob- 
tained from ([T]) by considering events £r involving pairs of components (Yj, Yh), j ^ h = 
1, . . . , g, of the random vector Y, i.e. 



n q—1 
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where fjh{-,-;9) denotes the marginal density of {Yj,Yh). The pairwise log likelihood 
is defined as pi{9) = log pL{9). The validity of using pairwise likelihoods to conduct 
inference about 6 can be assessed either from the theory of unbiased estimating functions 



(Godambe anc 


Kale 


Lindsay et al. 


2011) 



2005 



2011). The maximum pairwise likelihood estimate 9p is defined implicitly 



as the solution of the pairwise score equation 



9-1 



i=i j=i h=j+i 



dlogfjh{yij,yih? 
89 



0. 



Since Wi0[ps{9;Y)] = 0, the pairwise score function belongs to the class of unbiased 
estimating functions and 9p inherits the properties of M-estimators. Under regularity 
conditions assumed hereafter (see, e.g., Molenberghs and Verbeke, 2005), the maximum 



pairwise likelihood estimator is consistent and asymptotically normal, with covariance 
matrix given by the inverse of the Godambe information V{9) = H{9)~^J{9)H{9)~^, 
with J{9) = lEe[ps{9;Y)ps{9;Yy] and H{9) = lEe[-dps{9;Y)/d9'^]. 

Hypothesis testing and confidence regions for 9 can be obtained by using the anal- 
ogous of the Wald, the score and the log likelihood ratio tests. The pairwise likelihood 
counterparts of the Wald and score statistics are 



9yVi9)-\9p-9) and pWs{9) = psi9y Ji9)-^ps{9), (2) 



respectively, and both are asymptotically distributed as a chi-square random variable 
with p degrees of freedom. The pairwise log likelihood ratio 



pW{9) = 2 pi{9p) - pi{9) 



(3) 



converges in distribution to jyj=i^j{9)Zj, with Xj{9) eigenvalues of H{9) ^J{9) and 



Zj independent random variables having a standard normal distribution (Kent, 1982). 



The quantiles of the asymptotic distribution of ^ can be approximated by numerical 
algorithms (see, e.g., Imhof, 1961). The main drawback of tests and confidence sets 
derived from pW{9) lies in the fact that they might not enjoy the desirable large sample 
properties of their likelihood counterparts, as pW{9) is not asymptotically pivotal, i.e. 
its asymptotic distribution still depends on 9 through \j{9), j = 1, . . . ,p. 
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Approximate pivots can be obtained from pW{9) by suitable adjusting factors. A 
first statistic is obtained by a magnitude adjustment that forces the expected value of 
the asymptotic distribution of pW{6) to match the first moment of a chi-square random 
variable with p degrees of freedom (Geys et aL\ 1999). The resulting statistic is 



pW{9) 



(4) 



with Ki = Aj(6')/p, and its asymptotic distribution is only roughly chi-square as 

Ki corrects only the first moment of pW{6). Other moment-based adjustments can be 
considered. For instance, first and second moment matching gives the Satterthwaites 
adjustment (Satterthwaites, 1946) suggested in Varin (2008), whereas matching of mo- 



Further adjustments to pW{6) have been proposed by Chandler and Bate (2007): 



Wood 


(1989 


) and 


Lindsay et al. 


(2000 



pWM = pW{e)- 



and by Pace et al. (2011): 



pWinM=pW{e) 



-eyv{e)-\9p-e) 



ps{eyj{e)-^ps{e) 
ps{eyH{ey^ps{e)' 



(5) 



(6) 



The statistic ([s]) essentially stretches the pairwise log likelihood on the 6'-axis about 9p to 
ensure that the second Bartlett's identity holds. The statistic pWinv{d) can be derived 
from ([5]) by considering the formal relation (9p — 9) = H{6y^ps{6) + (9p(n^^). The main 
advantage of ([s]) and ^ over pWi{9) and other statistics derived from moment-based ad- 
justments is that they are asymptotically chi-square distributed and then asymptotically 
pivotal. 



1.2.2 Issues related to asymptotic variance estimation 

The matrices J {6) and H{6) determine both the convergence of statistics pWyj{6), pWs{0), 
pWi{6), pWcb{d), and pWinviO) to the central chi-square distribution and the quantiles of 
pW{9). In order to understand the way in which J {6) and H[9) affect the level error of 
tests and confidence sets derived from the aforementioned statistics, it is crucial to distin- 
guish two relevant scenarios in the pairwise likelihood framework. In the first one, pairwise 
likelihoods are used in place of the genuine likelihood function for computational conve- 
nience. Therefore a joint distribution for Y may be specified, and consequently either 
analytic expressions or Monte Carlo estimates for J {9) and H{9) can be worked out. In 
the second one, pairwise likelihoods are employed as approximations to the full likelihood 
function, i.e. only marginal bivariate distributions for sub- components of Y are speci- 
fied. In this case empirical counterparts of the elements of the Godambe information are 
needed. When dealing with independent observations, J{9) and H{9) can be consistently 
estimated by J{9) = n'^ Yl'I=iP-'^i^'^yi)P^i^'^yiV and H{9) = -n'^ Y^^^^dps{9;yi)/d9'^ , 
respectively. Otherwise, J{9) can be estimated by means of a window subsamphng esti- 



mator (Heagerty and Lele 19981 Heagerty and Lumley 2000), whereas the estimate of 



H(9) retains the structure of H(9). 
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The second scenario is far to be only a subtle distinction from the first one because 
to retrieve an accurate and stable estimate of J{6) and H{6) is still an open issue in the 



pairwise and, in general, composite likelihood framework (see, e.g, Varin et al. 2011, and 
references therein). In particular, large sample properties of pairwise likelihood statistics 
are affected both by the use of J (9) and H{0), and by the use of 9p in place of 9 in the 
computation of such estimates. When the sample size is moderate to small, replacing 
J{9) and H{9) with J{9p) and H[9p) slowdown the rate of convergence of the sampling 
distribution of pWw{9), pWs{9), pWi{9), pWcb{9), and pWinv{9) to the corresponding 
asymptotic references, leading to tests and confidence sets whose levels might be distant 
from the nominal ones. On the other hand, the estimates J{9) and H{9) improve, in 
general, the goodness of the approximation, therefore increase rejection and coverage 
accuracy of the derived tests and confidence sets, but need to be used carefully as H{9) 
might not be positive def i nite w hen considering values of 6' in a neighbour of 9p. As 
pointed out by Pace et al] (2011) replacing H{9) with H{9) is not an issue in hypothesis 
testing. However, it becomes relevant when considering non-null coverage probabilities 
of confidence sets, as statistics need to be evaluated at various values of 6* G 6. 

The above brief discussion reveals that the effect of estimating the elements of the 
Godambe information affects different properties of test and confidence sets derived from 
pairwise likelihood statistics. A rigorous mathematical treatment of such effects is difficult 
to assess, therefore in Section |3] a detailed account will be given through simulation 
studies. 



1.3 Some preliminaries on prepivoting 



Prepivoting has found important applications in reducing both the error level of tests and 
the coverage error of confidence regions. For the sake of simplicity a brief introduction 
to prepivoting is given by focusing on the former situation only. 

Consider the problem of testing the statistical hypothesis Hq : 9 = 9q, Hi : 9 9q, 9 E 
R. A statistical test at the level a, based on the pairwise log likelihood ratio, would reject 
Ho if pW{9y'' > qi-a, where pW{9)°'"' is the observed statistic, qi-a = Q~^{1 - a; Fe), 
and Q{x;Fg) = F[pW{9) < x]. 

In practice, the sampling null distribution of pW{9) is not known and the need of 
approximating qi-a leads to a test whose level is a only asymptotically. In finite samples 
the difference between the actual and the nominal level of the test mainly depends on 
the approximation of qi-a- Either asymptotic theory or nonparametric bootstrap can 
provide an approximation to the desired critical value. Nevertheless, as pW{9) is not 
pivotal the error level of the test would have the same order whatever approximation is 



adopted (Efron, 1982). 



When a non-pivotal statistic is considered, the bootstrap approach based on prepiv- 
oting can be used effectively to improve the asymptotic or the simple bootstrap ap- 
proximations of qi-a (Beran, 1987, 1988). Stemming from the test which rejects Hq 



if Q{pW{9)°^^; Fq) > 1 — a, the main idea of prepivoting is to move the attention from 
Q{x\ Fq) - which depends on 6' - to its null distribution function Qi{k] Fq) = P[Q{x] Fq) < 
k] which is uniform over the interval [0, 1]. Denoted with F some suitable estimate of Fg 



from which bootstrap samples are drawn, Beran (1987, 1988) shows that the bootstrap 
version Ql{k; F) = P*[Q*{x; F) < k] of the transformed statistic Qi{-; Fg) is less depen- 
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dent on 9 than Q*{x\ F) = P*[pW*{9) < x], where pW*{9) is the bootstrap version of the 
statistic and P* denotes probabihty with respect to F. 

Prepivoting can be iterated so that, at each iteration j, a bootstrap distribution 
Q*{u;F) = P*[Q*_i{k] F) < u] that is less dependent on 9, is built. In regular settings, 
it is possible to prove that, if Qi{-; Fq) is pivotal to order 0(n~*/^), then the distribution 
Q*{-;F) differs from the uniform random variable by an absolute error of magnitude 

0(n-*/2-i/2). 

The strength of prepivoting lies in its generality and in the opportunity to perform 
all the computations by Monte Carlo simulation. The generality of the method is paid at 
the price of a time consuming Monte Carlo simulation and, depending on the application 
area, the computational burden might relegate prepivoting to a theoretically attractive 
but practically unfeasible approach. In special cases, analytical prepivoting is possible 
(Beran, 1987, Section 3). Beran (1988) also discusses the possibility to reduce the compu- 
tational effort by using both analytical and mixed analytical-bootstrap approximations 
to prepivoting. However, in the present work these possibilities are not pursued since 
they would require the estimation of the elements of the Godambe information. 



2 Prepivoting in the pairwise likelihood framework 

2.1 The choice of the statistic and the resampUng plan 

In order to define a bootstrap test and confidence region there is the need to specify a 
suitable statistic and a sampling strategy that is consistent with the null hypothesis. 

Prepivoting statistics that are asymptotically pivotal already, as pWw{9), pWs{9), 
pWcb{9), and pWinv{9), would require a smaller number of bootstrap iterations to achieve 
a certain degree of accuracy than using a non-pivotal one. On the other hand, the the- 
oretical and computational advantage of bootstrapping asymptotically pivotal statistics 
would be annihilated by the collateral need of estimating the elements of the Godambe 
information. In fact, computing resampling-based estimates of J{9) and H{9) for each 
bootstrap sample would not necessarily cope with the issues related to variance estima- 



tion discussed in Section 1.2.2 Furthermore, prepivoting the aforementioned statistics 
would involve the computation of the maximum pairwise likelihood estimate. The models 
considered in the pairwise likelihood framework are usually rather complicated, thereby 
obtaining bootstrap versions of 9p could require an impressive amount of time. 

These considerations suggest that the use of a pivotal statistic is at odds with the need 
to obtain a resampling-based inferential procedure that is both accurate and reasonably 
fast. Instead, this trade-off may be avoided by focusing on a suitable non-pivotal statistic. 
In this paper it is proposed to use the unstudentized version of the pairwise score statistic 

pWU0)=n''ps{9ypsi9), (7) 

which converges to Aj(6')Z|, with Xj{9) the eigenvalues of J{9). As will be out- 

lined in the next section, the choice of ([T]) yields inferential procedures which achieve 
satisfactory levels of both accuracy and speed of computations. 

To compute the bootstrap null distribution of pWus{9) a suitable estimate F of Fg 
is needed in order to draw bootstrap samples y* = (y*, . . . ,?/*) consistent with the null 



6 



hypothesis. If the empirical distribution function F = was considered, y* would be 
reconstructed from y by using the uniform n-dimensional vector of resampling weights 
P 



n ^, . . . ,n ^). Nonetheless, this sampling plan may fail to supply the bootstrap 



null distribution of pWusifi) as it does not consider the possible invalidity of Hq (Hall 



and Wilson, 1991). To overcome this problem, it is possible to construct an estimate 
Fq centered at 6 from a vector of weights p{9) = {pi{0), . . . ,Pn{(^)) conceived to ensure 
that the bootstrap samples reflect Hq once 9 = 9q (Hall and Presnell, 1999). Here, it is 



suggested to obtain the functional form of the elements Pi{9) by minimising the forward 
Kullback-Leibler divergence between p and p{6) subject to J27=iPii^)P^i^Tyi) ~ 0- The 
analytic solution coincides with Owen's empirical likelihood formulation, therefore 



P^{0) 



1 



n{l + aOVps{9-y,)y 



where ^{9) e W solves ^'^^iPs{9;yi)/[n{l + ^{9yps{9;yi))] = 0. More details abo ut the 
derivation of ([sj) and the algo rithm used to obtain the root ^(9) can be found in Owen 
(flggol) andlHall and La Scalal (|T990|. 



The specific choice of Fg is primarily addressed by the need to obtain bootstrap 
samples that reflect the null hypothesis. In addition, it turns out that Fg enhances the 
effects of prepivoting by lightening the computational effort required in the Monte Carlo 



simulation, and will be clarifled in Section 2.4 



2.2 Computation of test and confidence set 

Let QUx; Fg) = P [pWUO) < x], Qusi{k] Fg) = P [Q„,(x; Fg) < k], and pW:!%9) be the 
observed value of pWus{9). The proposed a level test for Hq : 9 = 9q, Hi : 9 ^ 9q rejects 
the null hypothesis if 



pw:r{9) > {Q:,-Fg)-'[{Q:,,-Fg)-\i-a)], 

whereas the associated confldence set of level 1 — a for 6* is 

T^s = {9eQ: pWusi9) < (g:,; Fg)-'[iQ:,,; Fg)-\l - a)]} 



(9) 



(10) 



where QLls Fg) and <0«si('; Fg) are approximations to Ql^i-; Fg) and Ql^ii-; Fg), thereby 
bootstrap estimates of Q„s(-; Fg) and Qusi{-'i Fg). 

The estimates Qus{-] Fg) and Qusii''^ Fg) are usually obtained via Monte Carlo simu- 
lation: for the former, one outer level of 6 = 1, . . . ,B bootstrap replications is required, 
whereas the latter needs an additional inner level of m = 1, . . . , M bootstrap replications 
for each b. As the total number of computations equals B x M, some strategies have 



been proposed in order to lighten the computational effort. DiCiccio et al. (1992) propose 



to replace the inner level of bootstrap by using saddlepoint approximations to estimate 
Qusii,-] Fg). This approach is appealing but requires ad-hoc calculations for the model 
under consideration and is formally applicable in the smooth function of means model 



(Bhattacharya and Ghosh 1978). Lee and Young (1996) propose an algorithm embedding 



a stochastic stopping rule in order to reduce the number of inner bootstrap replications. 
Their algorithm is usually slightly less accurate when compared to the full-blown one 
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that entails B x M replications and requires the specification of some parameters to be 
ran. Nankervis (2005) suggests an approach which involves the use of a deterministic 
stopping rule in the inner level that allows to obtain the same results of the full-blown 
algorithm but with a smaller total number of bootstrap replications. 

In this paper the latter strategy is pursued and the resulting algorithm resembles 



Beran's original one (Beran, 1987) but it is modified to encompass the generation of 
samples according to the null hypothesis, i.e. F„ is replaced by Fg, and to embed a 
deterministic stopping rule. The main steps can be summarized as follows: 

0-Preliminaries: Let y = {yi, . . . ,yn) be the original sample and let a be the desired 
significance level for test ^ and confidence set (10). Evaluate and store the pair- 



wise score contributions ps{6;yi), i = 1, . . . ,n. Compute the weights Pi{9) of Fq 
according to (|8]), and attach to each element of the set of indices C = {1, . . . ,n} 
probability Pi{0); 

1- Outer level: For b = 1, . . . , B sample with replacement n elements from C obtaining 

the b-th new set of indices Cb. Compute the 6-th bootstrap version of pWus{0) given 
by pW*,.bi^) = n-\^i^^^ps{e; yi)V {^i^c^P^i^'^ V^)) ^^^^ store it along with A; 

Intermediate step: Sort the values pW*^.^^{6) into descending order, so that pW*^.^-^-^ 
and are the maximum and minimum values of the bootstrap repli- 

cates computed in the outer level, respectively. In an obvious notation C{p) is 
the set of indices associated to V^us-ib)^^)^ b = 1, . . . , B; 

2- Inner level I (full-blown): For the largest j = 1, . . . , + 1)J bootstrap values, 

use the corresponding score contributions indexed by in order to obtain a new 
vector of resampling weights p*{9) computed according to ([s]), and attach to each el- 
ement of probability i = 1, . . . ,n. For m = 1, . . . , M sample with replace- 
ment n indices from Cq) obtaining the m-th new set Cm- Compute pWus^Xi^) = 
-HE^ec^ Ps{0; y^)V{E^ec^ Vs{e- y,)) and A, = Zm=i M''! {pW:^sU0) < pW, 



n 



us;{j) 

Intermediate step: Sort A = (Ai, Aia{B+i)i) into descending order; 

3-Inner level II (stopping rule): For each of the remaining j = [Q;(i? + l)J +1, . . . ,B 

bootstrap values do the following: 

a) Use the corresponding score contributions in and obtain a new vector 

of resampling weights p*{6) and attach to each element of probability 
p*{e), i = l,...,n. 

b) Start the inner loop and at each iteration m* compute pW*g''l^,{6) and check 

whether 



^ m 

— J2 I {pW:iU(^) < +M-m*< A(L,(s+i)j) 

if this condition is satisfied stop further computations and go to Step a). 
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c) If all the M bootstrap iterations are carried out check whether 

1 ^' 

m=l 

if it is the case substitute A(^^a(B+i)\) with 7q), sort A into descending order 
and go to Step a). 

The outer and inner levels provide the desired bootstrap estimates of Qus{-]Fg) and 
iQusi] Fg)~^{l — a). In particular, 

1 ^ 
i=i 

and 

(QLi;^e)-'(i-«) = A(L,(5+i)j). 



2.3 Accuracy of test and confidence set 



In the following, the magnitude of errors entailed by ([oj) and (10) are provided under 
assumptions supplied in Hall (1992) and Lee and Young (2003). In order to ease the 
notation, the case p = 1 is considered as the results for generic p can be elicited with 
some minor modifications. 

Denote with k^, r = 2, 3, . . ., the cumulants of n~^/'^ps{0) (note K2 = J{0))- From the 



central limit theorem follows n ^/'^ps{6) ~ iV(0, K2) and n ^ps{9y 



1^2X1- 



Therefore the 



distribution function of pWusifi) 
Q 



n 



ps{6) may be expanded as 




+ 0{n- 



'IV. 



where g{-) 
and K3, 

respectively. In analogy, the bootstrap counterpart of (llll) is 



involves Hermite polynomials of order 3 and 5 which depend smoothly on 
$(■) and 0(-) are the standard normal distribution and density functions. 




X 
^2 




n 



-9 




X 
^2 



X 




K.2 



+ 0(n- 



(12) 



where §{■) has been obtained from (11) by replacing population cumulants with their 
bootstrap versions. Assuming that the difference between the estimates kr and their 
population counterparts is Op(n~^/^), the comparison of (11) and (12) shows that the 
bootstrap does not improve on the asymptotic approximation when considering a non- 
pivotal statistic. Moreover, Ql^g{-;Fo) is easily seen to be pivotal to order 0{n~^^'^). 
Finally, as the bootstrap is performed in a weighted fashion, i.e. samples are drawn 
according to Fg rather than to Fn, the following proposition can be derived by exploiting 



the results of Lee and Young (2003). 



Proposition 1. Under conditions in Lee and Young {2003} the difference between 
the actual and nominal levels of both ^ and (10) is 0(n~^^. 

The proof of Proposition 1 is omitted since it is sufficient to apply lemma (Al) and 



Proposition 2 of Lee and Young (2003) to (12). 
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2.4 Some remarks about the proposed approach 



Prepivoting the unstudentized version of the pairwise score statistic makes the use of 
nonparametric bootstrap appeahng in the pairwise hkehhood framework since the es- 
timation of the elements of the Godambe information is circumvented while obtaining 
fairly accurate inferential procedures that keep the computational burden under control. 
In the following, some key features of the proposed approach are briefly addressed and 
discussed. 

In first place, it is worth to provide an account of a slightly controversial point pursued 
in the former sections that spreads its implications in theoretical and practical aspects 
of the proposed prepivoting strategy. An Edgeworth view of the bootstrap state that 
resampling the non-pivotal statistic pWus{(^) rather than the asymptotic pivot pWs{6) 
would provide worse results in terms of accuracy of the derived tests and confidence sets. 
In fact, the relevance of asymptotics is dimmed because bootstrapping the latter statistic 
requires estimation of J{6) that is recognised to be an open issue in the composite like- 
lihood framework (see Section 1.2.2). Under these conditions, the bootstrap will likely 



be supplied with an unstable estimate of J{9) that may compromise the accuracy of 



the bootstrap (Hall et ai, 1989). As this effect can not be properly detected from the 



pertinent Edgeworth series, from a practical point of view it is considered more fruitful 
to violate the principle which favours to resample asymptotic pivots by relying on the 
bootstrap distribution pW*g{9), that automatically accounts for J (9), and by lowering the 
error level of the associated tests and confidence sets via bootstrap iteration. On the other 
hand, the lack of pivotalness of both pWus{9) and pW*g{9) implies that such statistics do 
not posses, in general, the desirable property of invariance under reparametrization, con- 
trasted to pW{9), pWs{9), a.ndpWinv{9). However, as pointed out by Pace et al. (2011), 



once that J{9) and H{9) are replaced by J{9p) and H(9p) the latter statistics lose exact 
invariance. 

In second place, the interdependence between computational costs and accuracy of 
test ([9]) and (10) need to be outlined. As a matter of facts, bootstrapping pWusi9) is the 
best choice for time saving. In order to form the bootstrap versions of pWus{9) at both 
the outer and inner levels, the bootstrap counterparts of the pairwise score contributions 
ps{9; Hi) are needed, z = 1, . . . , n. Since the value of 9 is fixed this allows to compute only 
once the pairwise score contributions and to reuse them by only sampling the indices in 
C and Cj. Furthermore, the use of pWus{9) avoids the computation of the maximum 
pairwise likelihood estimate B x M times, as would be required by using the statistics 
introduced in Section[l.2.1[ As a byproduct, the speed of the computations makes possible 



to choose the values of B and M not only on the basis of time constraints, but also to 
provide bootstrap estimates Qus{''^^e) and Qusii''^^e) that are reliable when considering 
critical values lying in the tail of the distribution of Qus{-', Fq). These considerations, that 
merely regard computational matters, must be further embedded into the resampling plan 



presented in Section |2.1| in order to better understand how accuracy of the proposed test 
and confidence set, claimed in Proposition 1, is achieved with a little computational 
expense. The results in Lee and Young (2003) state that if a statistic is pivotal to order 
0(n~*/^) and if p{9) is obtained by minimising the following divergence 



P(l-P) 



n 



Y.^npmY 



-00 < p < 00, 



i=l 



10 



then one bootstrap iteration sampling from Fq yields to a transformed statistic that is 
pivotal to order 0(r2~*/^~^) rather than 0(n~*/^^^/^) as would be obtained by sampling 
from Fn- (The suggested Fq is constructed from the vector p{9) obtained by minimising 
Dp subject to J27=iPii^)P^i^'^yi) ~ with p — )■ 0.) Therefore, since Qus{-',0) is pivotal 
up to 0(n~^/^), third order accuracy of ^ and (10) is reached with only one level of 
bootstrap iteration rather than two. 

In third place, performing weighted rather simple bootstrap has been shown to be 
both necessary in order to obtain the bootstrap null distribution of pWus{0) and to 
strengthen the effects of prepivoting. Some concerns may arise when computing the 
vector of resampling weights p{0), whose elements have the functional form provided by 



Owen (1988). In fact, depending on the sample size and on the dimension of 9, the 



convex hull condition might not be satisfied, resulting in a degenerate resampling vector 
p{6) which assign mass 1 to one unit (see, e.g., Owen, 2001). Experience from numerical 



investigations indicates that typically the occurrence of the convex hull issue is rather 
limited and can be regarded as a minor concern. 



3 Simulation studies 



3.1 Objectives 

In order to strengthen the soundness of the proposed approach, a simulation study has 
been conducted, serving two aims. The first aim is to provide an account of the accuracy 
of test ([9]) and of the associated confidence set (10) both in absolute terms and compared 
to the canonical and pairwise likelihood counterparts presented in Section 1.2. 1[ The sec- 
ond aim is to give numerical evidence of the motivations justifying this work by showing 
the impact of estimating the matrices J{d) and H{6) on null and non-null empirical rejec- 
tion and coverage probabilities of tests and confidence sets based on pairwise likelihood 
statistics. For this purpose the superscripts "n" and "e" will be used to denote statistics 
or quantiles computed by using the estimates J (9) and H{0), and those obtained by using 
J (Op) and H{Op), respectively. 

In the simulation setting the number of Monte Carlo trials have been set equal to 
20000, and the estimated quantile for test ^ and confidence region ( 10 ) has been obtained 
with B = M = 3000. 

The models from which data have been simulated, along with a summary of the 
associated results are described in the following section. 



3.2 Multivariate normal model 

As a first example, a rather simple multivariate normal model is considered. It serves the 
scope of comparing results of the application of the proposed approach with the u se of 
the considered competitors in a simplified setting where 9 = Op (Mardia et al., 2009) and 



J{6) and H{6) are available (Pace et al., 2011). 

The random vector Y is assumed to be distributed as a g-dimensional normal with 
mean (/z, . . . G R'' and compound symmetric covariance matrix E, having diagonal 
elements o"^ > and off-diagonal elements a^p, with p G (— l/(g — 1),1). The pairwise 
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log likelihood function for 6 = (/x, cr^, p) is 

pm = loga \og{l- p)~^^^^^—^SSw + 

q{q - 1)SSb + nq{q - - /i)^ 
2a2(l + p) 

where SSb = E»=i Yl=ii.yih - mf and S^h/ = ELill/i " ^i^h = XlLi Vih/^l and 
^ = Er=i Tl=iyihlnq. 

Samples of size n = 20 and dimension g = 10 have been drawn by setting the true 
parameter components /i = 0, o"^ = 1 and correlation coefficient p in {0.25,0.5,0.75}. 
For each sample pWus{0), w{0) have been computed as well as pW{6), tWw{Q), V^siO), 
pV^iiG)) V^cbiQ), and pWinv{0) and their counterparts using the estimates of J{0) and 



H{9) presented in Section 1.2.1 



Table [T] shows the empirical rejection probabilities for tests based on the aforemen- 
tioned statistics. The proposed test (|9| exhibits actual levels that are close both to the 
nominal ones and to those provided by the gold standard log likelihood ratio. Tests 
based on pairwise likelihood statistics computed by using the elements of the expected 
Godambe information or the estimates J{9) and H{9) exhibit levels close to the nominal 
ones, especially for pW{6) and pWinv{d). When J{6) and H{6) are replaced by J{Op) and 
H{9p) the error level of tests increases compared to the former situation, but for the ones 
based on pW{6) and pWi{6) that result to be more stable. 

In order to have a clue about the global reliability of test (|9| and tests based on 
pairwise likelihood statistics, it is useful to analyse the behaviour of the non-null empirical 
coverage probabilities of the associated confidence sets. In this analysis the parameter of 
interest is 6' = (cr^,p) and p is considered as known which allows an easier interpretation 
of the results by representing non-null coverage probabilities in contour plots. For each 
pair of (cr^,p) in an equally spaced 10 x 10 grid of points, probabilities are estimated via 
Monte Carlo simulation by generating samples with 6 = (a2 = l,p = 0.5),p = 0, n = 20, 
and q = 10. In Figure [T] are displayed the results for statistics reported in Table [T} The 
shape of contour plots provided by pWus{0), pWw{0), pWcb{0), and pWinv{0) are closer 
to that resulting from w{6) more than from those of pW{6) and pWi{d). Overall, all 
confidence sets exhibit null empirical coverages close to the nominal level. When pairwise 
likelihood statistics are computed by using J{9) and H{9) contour plots reveal that shapes 
become irregular and non-null empirical coverages do not decay to as moving away 
from (cr^ = l,p = 0.5), as would be expected, although null coverages for confidence sets 
derived from pW{6), pWi{6), and pWinv{G) remain quite close to the nominal level. The 
use of the estimates J{9p) and H{6p) provide contour plots whose shapes are similar to 
the ones obtained by using J {9) and H{9). However, in general, null empirical coverages 
are quite distant from the nominal level. 

3.3 Correlated binary data 

As a second example, a multivariate regression model with correlated binary response 
is considered. Besides being of practical interest than the previous one considered in 



Section 3.2 , the model presents more challenges because only the numerical evaluation of 



both the likelihood and pairwise likelihood function is possible. 
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Table 1: Multivariate normal model. Empirical rejection probabilities based on 20000 
Monte Carlo trials. Pairwise likelihood statistics denoted by the superscript "n" and "e" 
are computed respectively by using J {6) and H{9), and J {Op) and H{9p) 







p = 0.25 






p = 0.5 






p = 0.75 




a 


0.1 


0.05 


0.01 


0.1 


0.05 


0.01 


0.1 


0.05 


0.01 


w{e) 


0.100 


0.048 


0.009 


0.102 


0.050 


0.008 


0.103 


0.048 


0.009 




0.110 


0.052 


0.010 


0.111 


0.054 


0.009 


0.119 


0.060 


0.011 


pwm 


08Q 
0.108 
0.260 


042 
0.065 
0.193 


009 
0.022 
0.108 


111 

W . X X X 

0.136 
0.273 


062 
0.090 
0.207 


01 8 

W . W X (J 

0.043 
0.117 


1 75 

W. X 1 

0.205 
0.305 


1 22 
0.161 
0.241 


062 
0.100 
0.153 


nW (0) 

pw^{e) 

pwm 


0Q3 
0.186 
0.212 


051 
0.116 
0.151 


016 
0.039 
0.078 


092 
0.191 
0.204 


052 
0.120 
0.145 


016 
0.036 
0.080 


092 
0.193 
0.215 


052 
0.126 
0.159 


015 
0.036 
0.093 


r)W(f)) 

pWiO) 

pw\e) 


0Q8 
0.081 
0.080 


051 
0.039 
0.039 


01 
0.007 
0.007 


100 

w. xww 

0.082 
0.080 


050 
0.038 
0.038 


01 
0.007 
0.007 


095 
0.079 
0.077 


048 
0.039 
0.037 


01 
0.008 
0.008 


vWa (9) 

pwi^{e) 
pwi{e) 


0.105 
0.119 
0.147 


0.060 
0.065 
0.091 


0.017 
0.014 
0.034 


0.107 
0.109 
0.152 


0.060 
0.056 
0.095 


0.016 
0.010 
0.035 


0.114 
0.082 
0.154 


0.061 
0.038 
0.097 


0.015 
0.006 
0.030 


pw,t{e) 

pwm 
pwm 


0.091 
0.127 
0.240 


0.042 
0.065 
0.171 


0.008 
0.013 
0.086 


0.109 
0.152 
0.242 


0.056 
0.096 
0.168 


0.014 
0.040 
0.078 


0.161 
0.215 
0.259 


0.104 
0.163 
0.188 


0.046 
0.103 
0.096 


pW,nM 
pwuo) 
pwuo) 


0.091 
0.119 
0.238 


0.045 
0.055 
0.168 


0.010 
0.005 
0.082 


0.096 
0.082 
0.239 


0.046 
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0.062 
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Suppose that a binary outcome along with some relevant features are repeatedly 
measured on the same subject at q distinct temporal occasions. Let i/i be a g-dimensional 
vector and Xi he a q x p design matrix with ones in the first column which store the 
binary outcomes and the covariates for unit i, respectively, i = 1, . . . ,n. From a latent 
variable perspective (Renard et al, 2004), Yi = (Yn, . . . ,Yiq) can be thought of as the 



dichotomization of a continuous random vector Zj 



that is for some ^, 



Yij = 1 if Zij > ^ and Yij = otherwise. The random vector Z is assumed to be normally 
distributed with vector of means 7^ = Xj/3 depending on an unknown p-dimensional 
regression coefficient /3, and covariance matrix S, having diagonal elements > and 
off-diagonal elements a^p, with p G (— l/(g — 1), 1). 
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The log likelihood function for 9 = (/3,p) is 



J2^ogP{Y,=yf,9), (13) 



i=l 



where, for example, P(Fj = 1^; 6') = $q(cT ^'yi;p), with $g(-;p) the standard g-variate 
normal distribution function with correlation coefficient p and Iq a g-dimensional vector 



of ones. The evaluation of (13) becomes unfeasible as the number of observations over 
time increases, because $g(-; p) must be computed numerically. Resorting to the pairwise 
likelihood approach results in computational time saving as only bivariate integrals are 
involved. The pairwise log likelihood function for 6 is then 

n q—1 q 
1=1 j=l k=j+l 

where, as before, P(Yij = l,Yik = 1]9) = ^2{cr'^lij,<^'^lik; p), and and 7^^ are 
components of place j and of 7j, respectively. 

The simulation setting considers n = 15, g = 20, and p = 2 regression coefficients. 
After setting a = 1 and true parameter value to have components f3i = 0.5, (^2 = I 
and p in {0.25, 0.5, 0.75}, simulated data have been obtained accordingly to the following 
scheme. The design matrix Xj has been generated by considering q independent trials 
from a uniform random variable on the interval [—1, 1], whereas the binary outcome i/i 
has been obtained for unit i were obtained by drawing observations from Zi by setting 



,^ = (see, Renard et al, 2004, Section 3). 

For this model analytic expressions for J{6) and H{6) are not available and their 
corresponding empirical counterparts must be used to compute pairwise likelihood statis- 
tics. Table [2] shows the empirical rejection probabilities for test based on pWus{0), w{9) 
and pairwise likelihood statistics computed by using both J{0), H{6) and J{9p), H{9p). 
The actual levels of test ([9]) are close to the nominal ones and together with the full log 
likelihood ratio provides the best results. Also in this example, tests based on pairwise 
likelihood statistics exhibit a quite poorly behaviour, more marked when statistics are 
computed by using J{6p) and H{9p). Note that tests based on non-pivotal statistics 
pW{9) and pWi{d), whatever estimate of the elements of the Godambe information is 
used, outperform those based on asymptotically pivotal ones. 

More insights about the effects of estimating J {9) and H{6) on tests can be assessed 
by looking at the corresponding non-null empirical coverage probabilities of the corre- 
sponding confidence sets. The parameter of interest is 6' = {(32, p) and f3i is considered as 
known. Probabilities are estimated via Monte Carlo simulation by considering samples 
with n = 15, g = 20 and f3i = 0.5 for each pair (/32, p) in an equally spaced 10 x 10 grid of 
points. The true parameter value has components P2 = ^, P = 0.5, and nominal level is 
set to 0.95. In Figure [2] are displayed non-null empirical coverage probabilities for confi- 
dence sets obtained from statistics in Table 2 The use of estimates J{9) and H{9) leads 
to misbehaved contour plots that assign highest probability to values of 6 = {132, p) other 
than the true parameter value (see, in particular, the plot corresponding to pWinv{0))- 
On the other hand, resorting to J{Op) and H{6p) mitigates the problem without a remark- 
able deterioration of null coverages. However, whatever estimate is considered, none of 
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Table 2: Correlated binary data. Empirical rejection probabilities based on 20000 Monte 
Carlo trials. Pairwise likelihood statistics denoted by the superscript "n" and "e" are 
computed respectively by using J {6) and H{0), and J{9p) and H{6p) 
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these contour plots compare favourably with the one provided by w{6), neither in terms 
of shape nor in terms of null coverages. Confidence set Tus provides non-null empirical 
coverages that both are quite close to the ones of w{9) and outperforms uniformly those 
obtained from pairwise likelihood statistics. 

4 Final remarks 

Inferential procedures based on composite likelihood functions offer both flexibility in 
model specification and computational benefits. However, this potential is heavily com- 
promised by the need of estimating the matrices involved in the asymptotic variance of 
the maximum composite likelihood estimator. 

These problems are overcome by the fruitful application of resampling methods. 
Prepivoting the unstudentized version of the pairwise score statistic circumvent the esti- 
mation of the matrix J {6). Under suitable regularity conditions, the level of the derived 
test and confidence set has been shown to be third order accurate and the computational 
burden is kept under control. 

Simulation results in Section [3] confirm both that accuracy of tests and confidence 
sets derived from pairwise likelihood statistics deteriorates once that J{6) and H{6) are 
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estimated (see Section 1.2.2) and the benefits of bootstrapping a non-pivotal statistic (see 
Sections. 



First advocated by Aerts and Claeskens (1999), the opportunity of using bootstrap 
procedures in the composite hkehhood framework is far from being unexplored. However, 
the approach proposed in this work goes beyond the one of Aerts and Claeskens (1999) 
for general pseudo-log likelihood ratios. In first place, the bootstrap is performed in a 
nonparametric fashion, thus avoiding model assumptions that not always are affordable in 
the composite likelihood framework. In second place, Aerts and Claeskens (1999) propose 
to bootstrap pW{9) that, however, result to be less appealing from a computational point 
of view than pWus{0) as it requires the computation of the maximum composite likelihood 
estimate. Hence, the proposed approach opens an unexplored stream in the composite 
likelihood framework, where the use of nonparametric bootstrap leads to benefits in terms 
of both fiexibility and accuracy of the derived inferential procedures. 
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Figure 1: Multivariate normal model. Contour plots of non-null empirical coverage prob- 
abilities. Nominal level 0.95. From left to right: pairwise likelihood statistics computed 
respectively with J(^), E{Q); J {6), H{e)- J (Op), 11%). 
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Figure 2: Correlated binary data. Contour plots of non-null empirical coverage proba- 
bilities. Nominal level 0.95. From left to right: pairwise likelihood statistics computed 
respectively with J {6), H{e)] J {Op), H{ep). 
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